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Abstract. The results of simulations of extragalactic propagation of ultra-high energy cos¬ 
mic rays (UHECRs) have intrinsic uncertainties due to poorly known physical quantities and 
approximations used in the codes. We quantify the uncertainties in the simulated UHECR 
spectrum and composition due to different models of extragalactic background light (EBL), 
different photodisintegration setups, approximations concerning photopion production and 
the use of different simulation codes. We discuss the results for several representative source 
scenarios with proton, nitrogen or iron at injection. For this purpose we used SimProp and 
CRPropa, two publicly available codes for Monte Carlo simulations of UHECR propagation. 
CRPropa is a detailed and extensive simulation code, while SimProp aims to achieve accept¬ 
able results using a simpler code. We show that especially the choices for the EBL model 
and the photodisintegration setup can have a considerable impact on the simulated UHECR 
spectrum and composition. 
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1 Introduction 

The most energetic particles in the universe, the ultra-high energy cosmic rays (UHECRs), 
have been the subject of intense research for over fifty years. Although considerable progress 
in this field has been made in recent years with results from the two largest cosmic ray 
experiments, the Pierre Auger Observatory and the Telescope Array, the origin of these 
particles remains a mystery. The only aspect that has been widely agreed on is that the 
vast majority of the most energetic cosmic rays are charged particles (protons and/or other 
atomic nuclei) originating from outside our galaxy, and therefore they can interact with 
intergalactic background radiation and/or magnetic fields during their propagation to Earth, 
hence affecting their properties. 

At energies above 10 18 eV the cosmic ray spectrum presents some interesting features. 
One of these is the so-called “ankle”, at E ~ 5 x 10 18 eV, which is a flattening of the measured 
spectrum. Another feature is the flux suppression above ~ 2 x 10 19 eV [1, 2] . The flux suppres¬ 
sion may be a consequence of the well-known Greisen-Zatsepin-Kuzmin (GZK) effect [3, 4], 
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due to the interaction of UHECRs with the cosmic microwave background (CMB). Pure pro¬ 
ton models implementing the GZK effect, where the ankle is explained by proton interaction 
signatures as well, have been proposed extensively (see e.g. refs. [5, 6, 7]). However, the 
flux suppression could also be a signature of photodisintegration of nuclei by the CMB and 
the extragalactic background light (EBL), or of the maximum acceleration energy attainable 
by UHECR sources, as suggested in the context of the “disappointing model” [8]. The need 
for a vanishing proton component at E > 10 19 eV is driven by the increasingly heavy and 
pure mass composition as measured by the Pierre Auger Observatory [9]. In recent years 
combined investigations of the spectrum and composition measurements of Auger have shed 
a new light on the subject (see e.g. refs. [10, 11, 12, 13, 14, 15, 16, 17]). 

Possible models of UHECR sources can have many unknown parameters, such as their 
distribution, the maximum acceleration energy, the shape of the injection spectrum and the 
initial mass composition of cosmic rays. Therefore, to be able to thoroughly study UHECR 
source models and ascertain their compatibility with the available data, it is necessary to 
simulate the propagation of UHECRs in scenarios covering a large parameter space. To do 
this efficiently fast computational tools are required. There are a number of public codes for 
the propagation of UHECRs, including CRPropa [18, 19, 20], SimProp [21], TransportCR [22] 
and HERMES [23], as well as private codes (see e.g. refs. [24, 25, 26]), available for that 
purpose. 

The simulated UHECR spectra and mass compositions might depend strongly on poorly 
known quantities such as the spectrum and evolution of the EBL and photodisintegration 
cross sections, as well as on different computational treatments and approximations made 
in the different simulation codes. The main goal of the present work is to investigate how 
sensitive the simulations are to different models of the EBL, to different ways of treating 
photopion production, and to different photodisintegration implementations, as well as to 
compare the CRPropa and SimProp codes, aiming to understand how different computational 
treatments and approximations can affect the outcome of simulations. 

This paper is structured as follows: in section 2 the propagation of UHECRs in the 
universe is discussed; in section 3 the simulation codes are described; in section 4 the results 
of the comparisons are presented, for the cases of proton propagation and nuclei propagation 
separately; finally, in section 5 the results are discussed and the conclusions are presented. 

2 UHE cosmic ray propagation 

When propagating in the extragalactic space, UHECRs interact with photon backgrounds, 
namely the CMB and the EBL. The CMB has a blackbody spectrum with temperature at 
redshift 2 given by T(z) = To(l + z), with To = 2.725K. The EBL encompasses electro¬ 
magnetic radiation in a wide range of frequencies, from infrared to ultraviolet, and is poorly 
known mainly due to uncertainties concerning its time evolution 1 . 

The mean free path A for a particle with Lorentz factor T interacting at redshift 2 with 
a diffuse photon background of spectral number density n(e, z) for photons with energy e (in 
the laboratory frame) can be written as 
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x For details concerning the EBL models used in this work, refer to appendix B. 



where u(e') is the cross section for a given interaction between the cosmic ray and photons 
of energy e' = (1 — cos#)Te (in the nucleus’ rest frame), with 0 being the angle between 
the photon and the particle’s momentum in the laboratory frame. In the energy range of 
interest (E > 10 18 eV), the most important of such processes are photopion production, pair 
production and photodisintegration. 

Photopion production is the process in which a nucleon (free or bound to a nucleus) 
interacts with a background photon producing nucleons and pions (p+ 7 — > p+ 7T°, p+ 7 —> n+ 
7 r + , etc). Charged pions produce neutrinos and electrons ( 7 r + — > p + + 1 /^, p + —> e + + v e + v IJ ), 
whereas neutral pions produce gamma rays (7T° —> 7 + 7 ). The threshold energy for photopion 
production is ~ 3 x 10 19 (meV/e) eV, where e is the energy of the photon. The dominant 
photon background at ultra-high energies is the CMB (e ~ 0.7 meV), causing the so-called 
Greisen-Zatsepin-Kuzmin (GZK) effect, which induces a cutoff in the UHECR spectrum 
starting from energies around 3 x 10 19 eV. The GZK effect implies that almost all protons 
reaching Earth with energy greater than 10 2 ° eV must originate from within about 100 Mpc. 
For nuclei with mass number A, the threshold for photopion production is A times that for 
protons, but photodisintegration (see below) is possible at lower energies. 

Pair production is the process whereby electron-positron pairs are created due to the 
interaction of UHE nuclei with background photons (%X + 7 —> + e + + e~). It has 

a relatively short mean free path, but a very small fractional energy loss, thus being well 
approximated as a continuous energy loss (CEL) process. The threshold energy is ~ 5 x 
10 17 A(meV/e) eV, where A is the atomic mass of the nucleus. In the case of nuclei, the 
energy loss length (£ = T/(dT/dx)) scales as 7“^ = l~} otons Z 2 /A. 

Photodisintegration occurs when a nucleus is split into smaller parts due to interactions 
with photons [%X + 7 —> ^~ 1 X + n, + 7 —> z~\X + P- etc). Cross sections for this 
interaction are dominated by the giant dipole resonance (GDR) for photons with energies 
e' < 30MeV (in the nucleus rest frame). For 30MeV < e' < 150MeV quasi-deuteron (QD) 
processes occur, typically causing the emission of multiple nucleons. For e' > 150 MeV 
photodisintegration cross sections rapidly vanish and pion production takes over. Daughter- 
nuclei approximately conserve the same Lorentz factors from their corresponding progenitor 
nuclei as nuclear recoil is negligible, for the energy of interest is much less than the rest mass 
of the nucleus. 

Due to the expansion of the universe, all particles undergo adiabatic losses. In this case, 
the energy loss rate is given by 

- ^ = -H(z) = ^VOn + (l + z)*n A , (2.2) 

1 dx c c 

where Hq = 77(0) ~ 70km/s/Mpc is the Hubble constant at present time, D m « 0.3 is the 
density of matter (baryonic and dark matter), and Da ~ 0.7 is the dark energy density, in 
the standard cosmological model (ACDM). 

Another relevant process is nuclear decay. Unstable products of photodisintegration 
and photopion interactions can have lifetimes shorter than the typical propagation lengths 
involved, causing their decay along their trajectory to Earth. The most relevant nuclear 
decay processes for this energy range are a and decays, as well as nuclear dripping. 

A compilation of the energy loss length for different processes for 14 N and 56 Fe can be 
seen in figure 1 . 

Magnetic fields may affect cosmic ray observables measured on Earth, such as the spec¬ 
trum and mass composition. These effects may be due to the presence of intervening inter- 
galactic magnetic fields, the magnetization of the environment where the observer lies, the 
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Figure 1. Energy loss length for 14 N (left panel) and 56 Fe (right panel) at z = 0. Dotted-dashed 
lines correspond to interactions with the EBL, and dashed lines with the CMB. Green curves repre¬ 
sent electron pair production, purple photodisintegration, and orange photopion production. Double 
dotted-dashed lines are the energy loss length for adiabatic losses (Hubble radius). The thick black 
curve is the combined energy loss length, taking into account all processes. This plot was obtained 
using the interaction rate tables from CRPropa 3, with the EBL model of Gilmore et al. [27]. 


magnetization of sources and their surroundings, or any combination thereof. There are large 
uncertainties on the strength and distribution of magnetic fields, thus making it nontrivial to 
properly assess their impact on the UHECR spectrum and composition. On the other hand, 
in the case of sources uniformly distributed separated by a distance much smaller than the 
typical lengths involved (e.g. Larmor radius), it has been proven that the spectrum has a uni¬ 
versal shape regardless of the properties of magnetic fields [28]. A uniform source distribution 
is a reasonable assumption because the sources of UHECRs are not known, and considering 
the lower bounds on the density of sources estimated by the Pierre Auger Collaboration [29]. 
In refs. [30, 31, 32, 33, 34, 35] it was claimed that spectrum and composition may be affected 
by magnetic fields, although the rigidity at which this is relevant is not clear, ranging from 
10 15 V up to ~ 10 18 V. For E/Z > 10 18 V most of the aforementioned works predict small or 
negligible effects due to magnetic fields. 

3 Monte Carlo codes 

In the following the CRPropa and SimProp codes are briefly described. Attention is given to 
the considered models of the EBL and photodisintegration, as well as their implementation 
in the codes. 

3.1 The CRPropa propagation code 

CRPropa 3 [20] (see also refs. [36, 37, 38]) is a public 2 Monte Carlo code for propagating 
UHE nuclei, gamma rays and neutrinos in the universe. It includes all relevant interactions 
in the energy range of ~ 6 x 10 16 up to 10 22 eV, as well as many magnetic field environments 
and source distribution configurations. Three propagation modes are available, namely one¬ 
dimensional (ID), three-dimensional (3D), and four-dimensional (4D) modes. For the pur¬ 
poses of this work we will focus on the ID mode in order to compare with the SimProp code, 
which is limited to ID simulations. 

The energy loss processes and interactions implemented in this case are pair produc¬ 
tion, photopion production, photodisintegration, nuclear decay and adiabatic losses. Pair 

2 Code available at http://crpropa.desy.de. 
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production and adiabatic losses are treated as CEL processes, whereas photopion produc¬ 
tion, photodisintegration and nuclear decay are handled stochastically. The treatment of pair 
production follows ref. [39], and photopion production is handled by the SOPHIA code [40]. 

In CRPropa 3 photonuclear cross sections are obtained from the TALYS 1.6 code [41] 
using the parameters described in appendix A, in contrast to CRPropa 2 [19], which uses 
TALYS 1.0 with default parameters. TALYS has been used to predict photodisintegration 
products of all available exclusive channels: proton, neutron, deuterium, tritium, helium-3 
and helium-4 (alpha particle) and combinations thereof. Nuclei with A < 12 are not treated 
with TALYS. Instead, their photonuclear cross sections are compiled from various references. 
For 2 H, 3 H, 3 He, 4 He, 9 Be they are obtained from ref. [42], In the case of 3 H and 3 He they are 
scaled by factors 1.7 and 0.66, respectively, with respect to ref. [42], The parametrizations for 
9 Be is refitted to data, as shown in ref. [43]. Cross sections for 6 Li, 8 Li, 'Be, 11 Be, 8 B, 10 B, 
n B, 10 C and n C are taken from ref. [44], Cross sections for 'Li are obtained by interpolation 
of experimental data [45, 46]. For these nuclei, one proton is lost if Z > N, where N is the 
number of neutrons, one neutron if N > Z, or one of them with equal probability if N = Z. 
In CRPropa, a total of 185 isotopes (Z < 26, N < 30 with a lifetime r > 2s) and 2220 
disintegration channels are considered. Alternatively to the TALYS model, the total cross 
sections from ref. [44] can be used, while keeping the branching ratios from TALYS. 

Several EBL models are implemented in CRPropa 3, namely the ones by Kneiske et 
al. [47], Stecker et al. [48, 49], Franceschini et al. [50], Finke et al. [51], Dominguez et al. [52] 
and Gilmore et al. [27]. The implementation of photopion production, photodisintegration 
and electron pair production is based on tabulated mean free path data calculated with 
the comoving photon density n(e, 2 = 0) at redshift z = 0. By approximating a redshift 
independent spectral shape of the photon density, the following scaling relation for the mean 
free path A(r,z) at redshift z can be used (cf. [19]) 

H*) = (1 + ^) 3 {°oo n i £, nw e A((1 + ^ * = °) ( 3 -!) 

Jo n (e,0)de 

For photopion production CRPropa 3 alternatively provides redshift-tabulated values of the 
mean free path A(r, z) instead of the scaling relation. 

The treatment of nuclear decays is based on the NuDat 2.6 database 3 , which provides 
data for decay channels and nuclear lifetimes. In /3 + decays the absence of electron capture 
for fully ionized cosmic ray nuclei is accounted for in CRPropa. 

3.2 The SimProp propagation code 

The original version of SimProp, 4 described in ref. [21, 53] and henceforth referred to as Sim¬ 
Prop v2r0, can simulate the one-dimensional propagation of protons and nuclei in the absence 
of magnetic fields. All particles undergo adiabatic energy losses and pair production losses, 
treated as CEL with rates computed as in ref. [6]. Protons can lose energy via photopion 
production, approximated as a CEL with rates computed as in ref. [6]; pion production is 
neglected altogether in the propagation of nuclei. Nuclei undergo photodisintegration, which 
is treated stochastically according to the Puget-Stecker-Bredekamp (PSB) model [54] with 
Stecker-Salamon energy thresholds [55] (see also appendix A); all the ejected nucleons are 
treated as protons and the residual nucleus is treated as the stable isobar for the correspond¬ 
ing value of A. This version of the code computes interactions on the Stecker et al. EBL 

3 For details refer to the website http://www.imdc.bnl.gov/nudat2/. 

4 SimProp is available upon request to SimProp-dev@aquila.infn.it. 
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model [48, 49] (see also appendix B) as interpolated on a 2D grid of photon energies and 
redshifts, or a power-law approximation thereof [26, 56] 

An updated version, SimProp v2rl [57], treats pion production on the CMB as a stochas¬ 
tic process, both for protons and for nuclei, with rates computed using total cross sections 
from SOPHIA [40]; for nuclei, the rate is approximated as A times that for a proton with 
the same Lorentz factor; all photohadronic processes are treated as single-pion production, 
with branching ratio 1/3 for neutral pions and 2/3 for charged pions, in accord with isospin 
invariance. Also, the type of the nucleons ejected in photodisintegration is now randomly 
chosen and the residual nucleus is selected according to conservation of electric charge, so 
that photodisintegration and photopion production can now produce neutrons and unstable 
nuclei, though these are assumed to immediately undergo beta decay. Neutrinos produced 
in the decay of pions, neutrons and unstable nuclei are also tracked. 

SimProp v2r2 [58] also optionally considers pion production on the EBL, whose effect 
on UHECR fluxes is very small but relevant for cosmogenic neutrino production at energies 
below ~ 1 EeV. It also allows the user to choose the Kneiske et al. EBL model [47] using 
the photon energy at z = 0 scaled by a z-dependent factor as in CRPropa, as well as the 
ones used in SimProp v2rl. Moreover, it fixes a bug affecting the low energy tail of neutrino 
spectra, and uses a faster implementation to calculate interaction lengths resulting in much 
shorter computation times. 

In this work we use an extended version of SimProp (provisionally named SimProp v2r3) 
which will be released in the near future. SimProp v2r3 includes several new models of 
EBL, including Dominguez et al. [52] and Gilmore et al. [27] (see also appendix B), all 
interpolated on a 2D grid. Furthermore, it allows the user to specify different models of 
photodisintegration. In this work we use photodisintegration cross sections obtained from 
TALYS 1.6 with the settings described in appendix A, i.e. the same settings as used for 
CRPropa 3. 

We used the following approximation scheme to implement TALYS cross sections in 
SimProp v2r3, because while the uncertainties associated with the In A measurements [59] 
by the Auger observatory are in principle small enough to distinguish protons from helium- 
4, they are too large to distinguish consecutive intermediate nuclei, e.g. carbon-12 from 
carbon-13. Therefore it is important that UHECR propagation simulations accurately predict 
the number of protons and a-particles reaching Earth, but it is unnecessary to have the 
correct distribution of individual intermediate masses. In SimProp v2i'3, therefore, only two 
photodisintegration processes are implemented: nucleon ejection and a particle ejection; the 
interaction rates for these processes are taken to be the sum of all actual processes weighted by 
the number of nucleons and a-particles ejected, respectively. This ensures that the numbers 
of free nucleons and of a-particles at Earth, assuming that the interaction rates don’t change 
too rapidly with z or A, are reproduced in good approximation 5 * * * * * * . Since deuterium, tritium 
and helium-3 have very short disintegration rates, such ejectiles are treated as collections of 

5 For example, assume we have 14 N nuclei originating 70 Mpc away, and the only relevant process is 

14 N + 7 — > 12 C + p + n with interaction length 100 Mpc. A fraction exp(—0.7) « 50% of the nuclei will survive, 

and at Earth, for each 100 14 N nuclei injected, we will have in average 50 14 N nuclei, 50 12 C nuclei, and 100 

protons. If we chose to approximate this process as 14 N + 7 — > 13 C + p and 13 C + 7 — » 12 C + n with interaction 

length 50 Mpc each, a fraction exp(—0.7) 2 « 25% of the nuclei will survive, 2 exp(—0.7) (1 — exp(—0.7)) « 50% 

will interact once, and (1 — exp(—0.7)) 2 ~ 25% will interact twice, and at Earth, for each 100 14 N nuclei 

injected, we will have in average 25 14 N nuclei, 50 13 C nuclei, 25 12 C nuclei, and 100 protons. Both the 

number of protons and the average mass of the intermediate nuclei will then be well approximated, though 

the numbers of individual intermediate nuclides will be different. 
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free nucleons; also, since neutrons have a short decay length except at extremely high energy 
and even then the air showers they produce are indistinguishable from those of protons, we 
treat all nucleons as the same. Concretely, SimProp v2r3 considers these two processes: 


• nucleon ejection, with cross section 

(?N = ^ IT'N&n, 


' n 71 p 41 d 77 t 71 h 77 ck 


(tin)& tot i 


channels 


where un = n n + n p + 2nd + 3n t + 3nh, with the type of the ejected nucleon chosen 
at random with probabilities proportional to the proton and neutron numbers of the 
parent nucleus, and 

• alpha particle ejection, with cross section 



channels 


where cr nn n p n dnt n h n a is the exclusive cross section for ejecting n n neutrons, n p protons, ..., 
and n a a-particles as computed by TALYS. Finally, for faster computation of the interaction 
rates, &n and a a are fitted by a Gaussian for photon energies between the threshold energy 
and 30 MeV, and by a constant between 30 and 150 MeV. 

4 Comparisons 

In this section outputs of the simulation codes are compared. The comparisons are done for 
the case of pure composition at injection, with separate discussions for the case of protons 
(section 4.1) and nuclei (section 4.2) at the sources. For each case study, several scenarios 
are shown through comparisons of observables such as the flux at Earth and, in the case 
of nuclei injection, the average and variance of the logarithm of the mass. The observables 
are computed in log 10 (-E/eV) bins of width 0.1, from 17.5 to 20.5. The uncertainties in the 
choice of models are studied by computing (Ji~ Jj) / ((Ji + Jj) / 2 ) for the energy spectrum and 
(In A)i — (In A)j and u 2 (ln A)i — u 2 (ln A)j for the composition observables. These differences 
are presented together with the statistical uncertainties in the Auger data for the energy 
spectrum [60] and for the composition observables [59], the latter using the EPOS-LHC [61] 
hadronic interaction model. 

4.1 Propagation of protons 

In this section we compare the observed fluxes of protons, injected with a power law spectrum 
with spectral index 7 = 2.5 up to a maximum energy of 10 22 ' 5 eV, in order to assess the 
similarity between the simulation results for a wide range of initial energies. The source 
luminosity per unit of comoving volume is proportional to (1 + z ) 3 in the redshift range 
0 < z < 2.5. 

4.1.1 Stochasticity of pion production 

In order to study the effect of statistical fluctuations in pion production interactions, we 
compare SimProp results with the option -S 0 , where pion production is treated determinis¬ 
tically according to the CEL approximation, which has been frequently used in the literature 
(e.g. in ref. [ 6 ]), and with the option -S 1, where it is treated as a discrete interaction with 
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Figure 2. Comparison of observed proton fluxes simulated with pion production treated determinis¬ 
tically and stochastically. 


a stochastically sampled interaction point and energy loss. Here, the pion production is 
computed in the CMB only. 

The results are shown in figure 2. The two spectra are very similar, except that the 
peak at 10 19 ' 6 eV is broader in the stochastic simulation. The differences between the two 
fluxes are less than 10% at all energies and are only sizeable at E > 10 19 ' 5 eV, where existing 
measurements of the UHECR spectrum have large statistical uncertainties. 

4.1.2 Pion production on the EBL 

In order to study the effect of pion production on the extragalactic background light we 
compare the proton fluxes computed by SimProp taking into account the CMB and one of 
the following EBL models: Stecker et al. [48, 49], Kneiske et al. [47], Gilmore et al. [27], and 
Dominguez et al. [52], For a discussion and comparison of these models see appendix B. 

The resulting fluxes are shown in figure 3. Neglecting pion production on the EBL results 
in more than 10% higher cosmic ray fluxes at 10 19 ' 5 eV, but the difference is already visible at 
energies around 10 19 eV where the experimental data have small statistical uncertainties. The 
difference between EBL models is smaller, but is still at the level of the measured statistical 
uncertainties at energies around 10 19 ' 5 eV. 

Neutrino fluxes at Earth from the decay of pions and neutrons produced in the same 
scenario are shown in figure 3, both taking into account the CMB and the various EBL 
models, as well as the CMB only. It can be seen that even though the effect of EBL pion 
production on the proton flux is relatively small, it provides the majority of low energy 
(E < 10 1 ' 1 eV) cosmogenic neutrinos. 

4.1.3 Effect of different simulation codes: SimProp vs CRPropa 

We compare SimProp and CRPropa using the same EBL model (Gilmore et al. [52]), in order 
to investigate the effect of different propagation algorithms. In particular, the two codes 
use different step lengths for numerical integrations. In CRPropa the redshift dependence 
of the interaction rates for interactions with the EBL are approximated through a global 
scaling factor, whereas in SimProp rates are calculated at each redshift. In SimProp all 
photohadronic processes are approximated as single pion production isotropic in the center 
of mass frame, and with branching ratios from isospin invariance; neutron decay is treated 
as instantaneous. CRPropa treats photopion production of protons and neutrons differently, 
using the SOPHIA code to compute the energy dependent branching ratios and energy losses. 
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Figure 3. Comparison of observed proton fluxes (top) and cosmogenic neutrino fluxes (bottom) for 
different EBL models (Stecker, Kneiske, Dominguez and Gilmore), and for neglecting interactions on 


the EBL. 



proton injection, (1+z) 3 E~ 25 



Figure 4. Comparison between proton fluxes at Earth simulated using SimProp v2r3 and CRPropa 3. 


The simulated proton fluxes for both codes are shown in figure 4. The differences are 
small (< 10%), although they systematically depend on the energy, decreasing by about 
4% per energy decade. This results in a steeper spectrum for SimProp than for CRPropa, 
corresponding to a variation in the spectral index of about 0.02. 

4.2 Propagation of nuclei 

In this section we study the uncertainties related to the propagation of nuclei. The gen¬ 
eral scenario consists of identical sources emitting nuclei with a power law spectrum with 
rigidity-dependent cutoff dN/dE oc EP? exp(—E- m j/ZR cut ). The sources are homogeneously 
distributed in comoving volume, and the redshift range used here is 0 < z < 1. We con¬ 
sider two representative primary nuclides, nitrogen-14 and iron-56, and two representative 
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Figure 5. Effect of pion production on nuclei for soft nitrogen injection. 

injection spectra with associated rigidity-dependent cutoffs, 7 = 2 ,R cut = 10 20 V (“soft”) 
and 7 = l,i? C ut = 5 x 10 18 V (“hard”). (The soft injection is inspired by the best fit of the 
Auger spectrum and composition results reported in [10] in the case of pure iron injection. 
Hard injection of intermediate nuclei is found to be the best fit to Auger data from several 
authors [11, 13, 62] for an arbitrary mixed composition at the source.) We have simulated 
the resulting observables for all combinations of injection characteristics (hard iron, soft 
iron, hard nitrogen, soft nitrogen), for every propagation model studied. In the following, 
we show a selection of the investigated source characteristics and briefly describe qualitative 
differences for the cases that are not shown. 

4.2.1 Effect of pion production by nuclei 

In order to study the importance of pion production by nuclei, we compare fluxes and compo¬ 
sitions computed with SimProp for the provided options of neglecting and of considering pion 
production on bound nucleons. In the former (option -S -l), pion production for protons 
is approximated as a CEL process, whereas in the latter (option -S 1 ) it is treated stochas¬ 
tically. In both cases, the Gilmore et al. [27] EBL model and the PSB photodisintegration 
model are used. 

The results are shown in figure 5 for the case of soft nitrogen injection. It can be seen that 
the effect of pion production by nitrogen nuclei is negligible, except at energies E > 10 197 eV 
where the available measurements have large statistical uncertainties. This effect is even 
smaller in the case of iron primaries, since the energy threshold for pion production increases 
with nuclear mass, and in the case of hard injection, as the lower injection cutoff implies 
fewer primaries above the threshold. 

4.2.2 Effect of different EBL models 

To study the impact of different EBL models on the photodisintegration of nuclei, we have 
used SimProp with the two most up-to-date EBL models, namely Gilmore et al. and 
Dominguez et al.. The differences between these models are discussed in appendix B. In 
both cases the PSB photodisintegration model is used. The results are shown in figure 6 for 
hard iron injection and in figure 7 for soft iron injection. 
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Figure 6. Comparison of EBL models (Gilmore vs Dominguez) for hard iron injection. 








Figure 7. Comparison of EBL models (Gilmore vs Dominguez) for soft iron injection. 


The differences in the total energy spectra are significant, often exceeding 10%, resulting 
in softer spectra at Earth for the stronger Dominguez EBL model than for the Gilmore 
model. The differences in the spectrum are more clearly visible in the hard injection scenarios 
(more than 40% at E ~ 10 19 ' 3 eV), mainly due to different predicted numbers of low energy 
secondaries, while in soft injection scenarios low energy secondary protons are subdominant 
with respect to primary nuclei. Conversely, the differences in the average logarithmic mass 
are larger for hard injection, because the composition is dominated by secondary protons at 
low energy and primary nuclei at high energy with either model, whereas with soft injection 
the composition is more mixed and more model-dependent. 

The differences in partial spectra are mostly visible for low-energy intermediate mass 
secondaries of iron, because they are produced via repeated disintegration by the EBL. At 
high energies the partial spectra are in good agreement, because high-energy nuclei are mainly 
disintegrated by the CMB. The cases of nitrogen injection give very similar results for both 
hard and soft injection. 
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Figure 8. Comparison of PSB and TALYS photodisintegration models for hard nitrogen injection. 

4.2.3 Effect of photodisintegration cross sections I: PSB vs TALYS 

To investigate the effect of different photodisintegration cross section models, we compare 
SimProp simulations with both the PSB and the TALYS photodisintegration models. In 
both cases we used the Gilmore et al. EBL model. 

The results are shown in figure 8 for the case of hard nitrogen injection. The main 
difference between the PSB and TALYS models of photodisintegration cross sections is that 
only the latter includes channels where a-particles, rather than single nucleons, are ejected. 
The effect of these channels can be seen in the spectra at Earth, where the helium flux using 
TALYS largely exceeds the one obtained using PSB. In the case of soft nitrogen injection, 
differences in the total energy spectrum are about half of those for hard injection (because 
in that case secondaries are subdominant with respect to primaries even at low energies) 
and those in the average logarithmic mass are similar. The cases of iron injections result in 
very small differences (because channels ejecting a-particles are even more disfavoured with 
respect to those ejecting single nucleons than in the case of nitrogen). 

4.2.4 Effect of photodisintegration cross sections II: TALYS vs Kossov 

A third model for the photodisintegration cross sections is provided by the parametrization 
of Kossov [44], which is used in the GEANT4 code. Since the model does not parametrize 
partial cross sections, its total cross sections can be used in CRPropa in combination with 
the branching ratios from TALYS. We compare the results of CRPropa simulations using the 
TALYS and Kossov cross sections. In both cases the Dominguez et al. EBL model is used. 

The resulting fluxes are shown in figure 9 for hard injection of nitrogen and iron nuclei. 
Using the total cross sections from Kossov results in a higher level of photodisintegration 
for iron nuclei compared to TALYS, resulting in a difference in the spectrum of around 20% 
for E > 10 19 eV. For nitrogen the difference is much smaller in the energy range with low 
statistical uncertainty. The results are similar for both hard and soft injection scenarios. 
In our simulations the differences between TALYS and Kossov are small, compared to that 
between TALYS and PSB, or between using different EBL models. This is due to similar 
total cross sections and because the same branching ratios (from TALYS) are used in both 
cases. 
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Figure 9. Comparison of TALYS and Kossov photodisintegration cross sections for hard injection of 
nitrogen (left) and iron (right) nuclei. 

4.2.5 Effect of photodisintegration cross sections III: TALYS with rescaled a a 

In this section we investigate the influence of photodisintegration channels where a-particles 
are ejected on the spectrum and composition observables. To this end, we ran SimProp 
simulations using TALYS photodisintegration where all values of a a (defined in section 3.2) 
were scaled by a factor a a = 1.0 (unsealed), 0.3, 0.1, and 0.0 (a ejection disabled). In each 
case the Gilmore EBL model was used. 

The reason for this kind of analysis is the lack of cross section measurements for the 
a-particle ejection: in the data sets used in the TALYS code the only measurements for this 
photodisintegration channel (for nuclei with 10 < A < 56 and photons below 30 MeV in the 
nucleus rest frame) are the ones of 12 C ( 12 C + 7 —> 3a, 12 C + 7 —>• p + a + ' Li) and 16 0 
( 16 0 + 7 —>• 4a). Moreover, TALYS seems to overpredict these measurements, as shown in 
appendix A. It is thus worthwhile to understand what the impact of these uncertainties is 
on the observables. Results are shown in figure 10 and following. 

The effects of changing this poorly known quantity on the energy spectrum in the case 
of hard nitrogen injection can be very large, over 50%, resulting in softer spectra at Earth 
the larger a a is. In the case of soft nitrogen injection these differences are about half as large, 
as there are more primaries at low energy and more secondaries at high energy than for hard 
injection. Conversely, the effects of this scaling on the average logarithmic mass in the case 
of soft nitrogen injection are large at all energies, whereas those for hard nitrogen injection 
are similar at low energies but negligible at high energies, where the composition is almost 
purely primary nitrogen regardless of the value of a a . In the case of iron injection, the effects 
of the scaling are smaller, because for heavy nuclei a-particle ejection is strongly disfavored. 

4.2.6 Effect of different propagation codes: SimProp vs CRPropa 

In order to study the effect of the different simulation algorithms on the propagation of nuclei, 
we compare fluxes and composition observables computed by SimProp and CRPropa, both 
using the Gilmore EBL model and TALYS photodisintegration cross sections. 

The results for hard nitrogen injection are shown in figure 13. It can be seen that once 
the two codes are used with the same models for the EBL spectrum and photodisintegration 
cross sections, the remaining differences due to the different approximations used in the 
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Figure 10. Effect of scaling of TALYS cross sections for a-particle ejection for soft nitrogen injection. 
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Figure 11. Effect of scaling of TALYS cross sections for a-particle ejection for hard nitrogen injection. 


hard iron injection 


1e+38 



a„=1.0 — 
a a =0.3 — 
1e+37 £ a a =0.1 - 
a„=0. 


1e+36 r 


18.5 19 19.5 20 20.5 

log 10 (E/eV) 



log 10 (E/eV) 



|> 0.5 

1“ 0 
< 

3 -0.5 

-1 

17.5 18 18.5 19 19.5 20 20.5 

log 10 (E/eV) 



Auger stat. unc. 
a a =T.O v| a a =0.3 
a a =1.0 vsa„=0.1 
a a =1.0 vs a a =0.0 

mm 


- 


17.5 18 18.5 19 19.5 20 20.5 

log 10 (E/eV) 


0.6 


1 1 B 

- 1-1 - 1 - 

0.4 



“ 

0.2 







$. 

0.2 


'% / 

Auger stat. unc. ■■ 


•BWy* 

a„=T.O vs a„=0.3 --- ~ 

0.4 


a a =1.0 vs a a =0.1 — _ 
a^=1.0 vsa a =0.0 

0.6 


-1-1- 


17.5 18 18.5 19 19.5 20 20.5 

log 10 (E/eV) 


Figure 12. Effect of scaling of TALYS cross sections for a-particle ejection for hard iron injection. 
















































































hard nitrogen injection 








Figure 13. Effect of different propagation codes for hard nitrogen injection 


algorithms are small (of the order of 10% or less except at the highest energies), although 
larger than the statistical uncertainties on the energy spectrum. The cases of soft nitrogen 
injection and iron injections result in similar or smaller differences. 

Possible reasons for the remaining differences include the simplified treatment of pho¬ 
todisintegration in SimProp described in section 3.2, the simplified redshift-dependence of 
photodisintegration on the EBL in CRPropa, described in section 3.1, or the different cross 
sections for light nuclei (cross sections as listed in section 3.1 for A < 12 in CRPropa, PSB 
cross sections for A < 5 in SimProp ). 

5 Discussion 

We studied the propagation of protons and its effect on energy spectra at Earth with fixed 
assumptions for the injection spectrum at the source. An analysis concerning the propaga¬ 
tion of nuclei has been done for two representative primary nuclides and two representative 
injection spectra. One of the main purposes of this study was to quantitatively estimate the 
impact of different propagation models in the observables. The observed differences in the 
spectrum and composition at Earth could imply variations in the source models predicted 
by several authors such as the ones proposed in refs. [11, 13, 62]. The differences in the 
calculated observables are shown together with the Auger statistical uncertainties in order 
to assess whether they can possibly have sizeable effects on the results of fits to Auger data. 

The effect of different implementations of photopion production has been tested. We 
have quantified how the approximation of photopion production as a continuous energy loss 
process, as done e.g. in refs. [6, 7], affects the measured spectrum and mass composition, 
compared to the stochastic treatment. There are visible differences at energies > 10 19 ' 5 eV, 
though they are below 10%. These differences tend to decrease for heavier nuclei, due to the 
dependence of the energy threshold for photopion production on the mass of the particle. 

Differences in the spectral density and redshift evolution of the EBL models are expected 
to have an impact on the studied observables. Most of the latest EBL models agree on 
the overall shape of the spectral density, particularly at z = 0, and in the ultraviolet and 
visible regions. Nevertheless, there are some differences that impact the propagation of 
UHECRs, which is mainly affected by the infrared region, affecting both the spectrum and 
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mass composition. These small differences tend to increase with redshift, thus impacting 
the study of cosmogenic neutrino fluxes. As a consequence, one needs to constrain the 
uncertainties due to the choice of EBL models by using, for example, multiple models or the 
available upper and lower bounds (see ref. [63] for a review). 

In the case of proton injection, among the EBL models considered, the largest bump in 
the UHE proton spectrum at Earth, apart from the one due to the CMB, is obtained using 
the Stecker et al. EBL model [49] whose intensity in the infrared peak is lower compared to 
the other models studied (see appendix B, fig. 15). The choice of the EBL model does not 
significantly affect the spectrum measured at Earth in this particular case, since the CMB 
dominates over the EBL. However, it is extremely important for studies involving cosmogenic 
neutrinos at PeV energies, since it is directly involved in their production. Nevertheless, one 
should bear in mind that although the flux of neutrinos produced due to proton interactions 
with photon fields can be affected by this choice, this uncertainty is smaller than the uncer¬ 
tainties in the cosmological evolution of sources [64]. In the case of nuclei, due to the giant 
dipole resonance, at energies ~ lOEeV photodisintegration can occur via interaction with 
photons with energies ~ 10 — 100 meV in the laboratory frame. For this reason, the EBL 
plays a fundamental role in the propagation of nuclei. In the particular case of the Gilmore 
et al. [27] and Dominguez et al. [52] EBL models the differences in the energy distribution 
of photons can result in observed differences of up to 40% in the spectrum, and 30% in the 
average logarithm of the mass number, exceeding statistical uncertainties from Auger, thus 
being non-negligible. As expected, differences are visible up to ~ 10 20 ' 2 eV in the case of iron 
injection, while in the case of nitrogen injection they are visible up to 10 19 ' 7 eV, due to the 
lower value of the threshold for the photodisintegration in the CMB. 

In ref. [25] the authors compare three EBL models, namely the ones from refs. [65, 66, 
67], and claim that there are considerable differences among them only for particularly heavy 
nuclei. Their conclusion is based on the comparison of energy loss lengths for photodisinte¬ 
grations assuming these EBL models at a fixed redshift (z = 0). Moreover, the EBL models 
compared in ref. [25] appear to have a larger impact at the lowest energies, where indeed the 
energy loss length of photodisintegration in the EBL, which at these energies is mainly due 
to the UV peak of the EBL density distribution, is much larger than the one for adiabatic 
losses. The IR peak of its distribution is then the main contributor to photodisintegration in 
the EBL. As shown in fig. 15, the discrepancies among different models are larger for larger 
wavelengths, making a discussion of their influence on the observables necessary. Further¬ 
more, the spectral energy distribution of the EBL evolves with redshift, making the cosmic 
ray spectrum a better observable to study the effects of EBL models than the energy loss 
length only. 

Different approaches in the propagation algorithms, as implemented in CRPropa and 
SimProp, have been investigated. In the case of protons the differences do exceed the sta¬ 
tistical uncertainties of Auger in the energy range between 10 18 ' 5 eV and 10 195 eV, although 
they are always smaller than 10%. The different approximations adopted for the treatment 
of photopion production are likely to be the main reason for the discrepancies. We have also 
compared the two simulation codes for the case of nuclei, adopting the same EBL model 
and photodisintegration cross sections. The differences are of the order of 10%, except for 
energies above 10 2 ° eV, in which case this number can reach 15%. Possible reasons for these 
differences are the simplified treatment of the photodisintegration in SimProp , the simplified 
treatment of the EBL redshift evolution in CRPropa, or the different cross sections used for 
light nuclei, besides the differences in the treatment of protons, which may be produced via 
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photonuclear interactions in the case of primary nuclei. 

A substantial part of the present work is dedicated to the discussion of the impact of 
photodisintegration cross sections on observables. The default version of SimProp uses the 
PSB photodisintegration model, which does not include channels whereby a particles are 
ejected, as TALYS does. The consequences of this simplification can be easily seen in the 
measured flux of individual species: there is a significant enhancement in the flux of helium 
nuclei and a reduction in the flux of intermediate mass nuclei such as nitrogen. This results 
in a decrease of (In A) in the energy range between ~ 10 18 ' 5 eV and 10 19 ' 5 eV. The lack of 
measurements of photonuclear cross sections in the data sets used by TALYS, and the fact 
that TALYS overpredicts some of the measured cross sections, motivated the study of the 
effect of rescaling the a production rates. Moreover, the uncertainties associated with (In A) 
by Auger are, in principle, small enough to allow us to distinguish between protons and 
helium nuclei. Therefore, it is important that UHECR propagation simulations accurately 
predict the number of protons and a-particles reaching Earth. The ejection of a particles is 
less frequent in the case of iron than in the case of nitrogen, and therefore the effects of this 
scaling in the iron injection scenarios are smaller than in the nitrogen cases. 

The parametrization of cross sections used in GEANT 4 [44] together with the branching 
ratios from TALYS is compared with the TALYS cross sections in CRPropa. Differences are 
small compared to those arising from using different EBL models, and from including the 
a-particle ejection. 

As a general remark, the differences in the spectra are often more visible in the hard 
injection scenario (7 = l,R cut = 5 x 10 18 V) than in the soft one (7 = 2,R cut = 10 2 ° V), 
because they are mainly driven by the ratio of primary to (low energy) secondaries, which 
is smaller in the soft than in the hard injection scenarios. On the other hand, considering 
the composition observables only ((In A) and cr(lnA)), the differences are smaller for harder 
injection spectra. This can be interpreted in terms of the low rigidity cutoff in the hard 
injection scenario, which implies a virtually negligible contribution of lighter nuclei at the 
highest energies, allowing only the residual primaries of heavier nuclei to arrive at Earth. 

In refs. [11, 13, 62] it has been argued that hard spectral indices (7 < 1.6) and lower 
maximal energies (E max ~ Z x 5 x 10 18 eV) are needed in order to simultaneously fit the 
spectrum and composition measured by Auger. Assuming a rigidity dependent cutoff for the 
injection spectrum, the cutoff at such low energies arises from the need to have a vanishing 
light component at the highest energies. As a consequence, the hard spectral indices are 
required to obtain the right ratio between residual primaries and lighter masses at Earth. 
In fact, as we have shown in the plots of the spectra, for sources injecting nuclei with a 
fixed maximum energy E max , the harder the injected spectrum is, the larger are the fluxes of 
secondary nuclei with respect to the residual primaries. Therefore, a detailed understanding 
of photonuclear processes is important in order to address this issue, since the flux of different 
particle species arriving at Earth may be significantly under- or overestimated when adopting 
different approaches. 

Understanding the impact of EBL models is also important. In fact, we have demon¬ 
strated that the enhancement of the efficiency of photodisintegration, as it occurs when the 
a-channel is included, has qualitatively the same effect as considering the most intense EBL 
model. One can then argue that a scenario in which the most intense EBL model is cho¬ 
sen together with a photodisintegration model that includes the a-channel will result in an 
extremely efficient production of secondary nuclei. Taylor et al. [16] assume that the un¬ 
certainties in the EBL distribution and in cross sections relevant for nuclei propagation are 
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unlikely to qualitatively impact the conclusions on the fit of the Auger data. Nonetheless, 
we have quantitatively estimated the influence of some combinations of EBL and photodis¬ 
integration models on the observables at Earth, showing that even with the simplified choice 
of pure injection the differences often exceed the statistical uncertainties of Auger. 

6 Conclusions 

In the present work we have discussed how sensitive the UHECR propagation is to un¬ 
certainties in the EBL spectrum and in the photodisintegration models. This has been 
investigated using two different simulation codes, CRPropa and SimProp , which have been 
directly compared in order to understand the effect of different computational treatments in 
the observables at Earth. Our results suggest that uncertainties in the scaling of a-channels 
related to the ejection of a-particles is the dominant source of uncertainties amongst all 
studied parameters, including different EBL models and photodisintegration cross sections. 
Also, the energy spectrum at Earth is more sensitive to the uncertainties in propagation 
in scenarios with hard injection spectra, whereas the measured mass composition is more 
model-dependent for soft injection spectra. 

To summarize, we have shown that different choices of parameters such as photonuclear 
cross sections, EBL model and computational treatment, can have a considerable impact in 
UHECR observables such as the spectrum as composition. The present work could be used 
as a quantitative estimation of uncertainties due to propagation in such interpretations. 
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A Photodisintegration cross sections 

One of the most important processes in the propagation of ultra-high energy nuclei through 
diffuse background radiation is photodisintegration, where a nucleus absorbs a photon and 
then ejects one or more fragments (most commonly single nucleons, but also a-particles 
or multiple nucleons). The experimental data about the cross sections of such processes 
are limited: for many nuclei only measurements of the total photoabsorption cross section 
and/or of single neutron ejection are available, mainly due to the difficulty in detecting 
outgoing charged particles. For this reason phenomenological models are needed in order to 
implement these processes in UHECR propagation simulation codes. 
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Isotope 

E 0 [MeV] 

C7 0 [mb] 

To [MeV] 

E 1 [MeV] 

o\ [mb] 

Ti [MeV] 

Source 

C-12 

22.70 

21.36 

6.00 




Atlas 

N-14 

22.50 

27.00 

7.00 




Atlas 

0-16 

22.35 

30.91 

6.00 




Atlas 

Na-23 

23.00 

15.00 

16.00 




Atlas 

Mg-24 

20.80 

41.60 

9.00 




Atlas 

Al-27 

21.10 

12.50 

6.10 

29.50 

6.70 

8.70 

RIPL-2 

Si-28 

20.24 

58.73 

5.00 




Atlas 

Ar-40 

20.90 

50.00 

10.00 




Atlas 

Ca-40 

19.77 

97.06 

5.00 




Atlas 

V-51 

17.93 

53.30 

3.62 

20.95 

40.70 

7.15 

RIPL-2 

Mn-55 

16.82 

51.40 

4.33 

20.09 

45.20 

4.09 

RIPL-2 


Table 1 . Giant dipole resonance parameters used with TALYS (as parameters for the Kopecky-Uhl 
generalized Lorentzian model of the El-strength function): peak energy Ei, peak cross section Oi 
and width Ej for resonances with a single (i = 0) or a split peak (i = 0,1). Default values from the 
RIPL-2 database are replaced, if available, with the total cross section parameters from the atlas of 
GDR parameters. Note that for isotopes not listed, as well as for higher order contributions, TALYS 
uses a compilation of formulas listed in [72]. 


The model proposed by Puget, Stecker and Bredekamp [54] includes a restricted list 
of nuclides (with one isobar for each A from 2 to 4 and from 9 to 56), and approximates 
the cross sections for one- and two-nucleon ejection for photon energies in the nucleus rest 
frame 2 MeV < e' < 30 MeV as Gaussians, and cross sections for multi-nucleon ejection 
for 30 MeV < e' < 150 MeV as constants, with tabulated branching ratios for the possible 
number of ejected nucleons. The exception is beryllium-9, for which the only photodisinte¬ 
gration channel is into two nuclei of helium-4 and one proton. A refinement of this model by 
Stecker and Salamon [55] uses the kinematic threshold for each process instead of 2 MeV as 
the lower limit of integration. Throughout this work, by PSB model we refer to this refine¬ 
ment. The PSB model makes no distinction between ejected protons and neutrons; when it 
is used in SimProp, the corresponding branching ratios are taken to be proportional to the 
number of protons and neutrons in the parent nucleus. Also, channels involving the ejection 
of fragments other than single nucleons (e.g. a-particles) are neglected. 

A more sophisticated model is provided by the nuclear reaction program TALYS [41]. 
It allows to compute cross sections for all exclusive photodisintegration channels, describing 
the ejection of protons, neutrons, deuterons, tritons, helium-3 and helium-4 nuclei, and any 
combinations thereof. A preliminary version of TALYS was used by Khan et al. [68] for an 
exhaustive comparison to the available experimental data. In their comparison TALYS was 
used with the giant dipole resonance parameters compiled in the atlas of GDR parameters 
in ref. [69]. In contrast, the publicly available versions of TALYS take these parameters from 
the RIPL-2 database [70] by default. In this work we make use of the former parameters, if 
available [71] (see table 1 for a complete list), as the resulting cross sections (listed as “TALYS- 
1.6 (restored)” in figure 14) are in much better agreement with the available measurements. 

Yet another model is that by Kossov [44] used in the Geant4 software. It describes the 
total cross sections for photodisintegration, as well as at higher energies pion production. 
Since branching ratios of individual disintegration channels are not modeled, these are taken 
from TALYS, when the model is used in CRPropa. 
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Figure 14. Photodisintegration cross sections for total absorption by silicon-28 (left) and a-particle 
ejection from carbon-12 (right) as predicted by various models. The measured data (yellow circles) are 
from ref. [73] and [74] respectively. TALYS (default) refers to using default TALYS settings, TALYS 
(restored) to using the GDR parameters listed in table 1. 


In figure 14 we compare cross sections predicted by the three models we used in this 
work with available measured data for the total photoabsorption cross section of silicon- 
28 and for the ejection of an a-particle from carbon-12 (with subsequent near-immediate 
decay of the residual beryllium-8 into two more a-particles). For completeness, we also 
show cross sections computed by the publicly released versions of TALYS using their default 
settings; in particular, TALYS-1.0 with default settings is the photodisintegration model used 
in CRPropa 2. It can be seen that TALYS with the parameters used in the original paper 
most closely reproduces the total cross section data, but PSB and Kossov also give acceptable 
results, whereas TALYS used with its default settings predicts much broader and lower peaks 
than observed. On the other hand, all versions of TALYS (especially TALYS-1.0) largely 
overpredict cross sections for a-particle ejection, as does the Kossov model with TALYS 
branching ratios (though at higher photon energies), whereas the PSB model neglects it 
altogether. Note that for kinematical reasons the cross sections at the lowest photon energies 
are the ones most relevant for UHECR propagation. 

B Models for extragalactic background spectrum 

The spectrum of the diffuse extragalactic background radiation spans over 20 decades in 
energy, from radio waves up to the high-energy gamma ray photons. It consists of light 
emitted at all epochs, modified by redshifting and dilution due to the expansion of the 
universe. The cosmic microwave background (CMB), the relic blackbody radiation from 
the Big Bang, is the dominant background field, followed by ultraviolet/optical and infrared 
backgrounds (extragalactic background light, EBL). In this work, several models for EBL 
have been used; these models are included in the simulation codes for propagation with 
different choices for considering how their spectral energy distribution evolves in redshift. 

The understanding of the spectral energy distribution and redshift evolution of the EBL 
requires studying the sources responsible for its production. Several different techniques are 
used for this purpose. Kneiske et al. [47] report the present-day background intensity using 
detailed information from galaxy surveys about global quantities as the cosmic star formation 
rate. The work by Dominguez et al. [52] is observationally based on multiwavelength data. 
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Figure 15. Intensity of the EBL at z = 0 (left) and z = 1 (right). Wavelengths are given in local 
physical lengths at the given redshift, and densities are per unit of comoving volume. 


Other authors (as for example Stecker et al. [48, 49]) use “backward evolution” of the present 
day galaxy emissivity. On the contrary, “forward evolution”, which begins with cosmological 
initial conditions and follows a forward evolution with time by means of semi-analytical 
models of galaxy formation, is used in Gilmore et al. [27]. 

In figure 15 the intensity of the EBL at z = 0 (left) and z = 1 (right) as a function of 
wavelength is shown, as predicted by the models used in this work (Gilmore 2012 [27] and 
Dominguez 2011 [52]) and by the default EBL models used in SimProp v2r0 and CRPropa 2 
(Stecker 2005 [48, 49] and Kneiske 2004 [47] respectively), as well as the Franceschini 2008 
model [50] for comparison. It can be seen that all recent EBL models are in good agreement 
concerning the EBL spectrum in the UV and optical region in the local universe, but they 
still largely differ in the far IR region and at high redshifts. (Note that due to the 1/e 2 
factor in eq. (2.1), the far IR region is the most relevant to UHECR propagation, and the 
UV region has very little impact even for the Stecker 2005 model at high redshift where it 
largely exceeds other models.) 

References 

[1] Pierre Auger Collaboration, J. Abraham et al., Observation of the suppression of the flux of 
cosmic rays above 4 x 10 19 eV, Phys.Rev.Lett. 101 (2008) 061101, [arXiv: 0806.4302]. 

[2] HiRes Collaboration, R. Abbasi et al., First observation of the Greisen-Zatsepin-Kuzmin 
suppression , Phys.Rev.Lett. 100 (2008) 101101, [astro-ph/0703099], 

[3] K. Greisen, End to the cosmic ray spectrum?, Phys.Rev.Lett. 16 (1966) 748- 750. 

[4] G. Zatsepin and V. Kuzmin, Upper limit of the spectrum of cosmic rays , JETP Lett. 4 (1966) 
78-80. 

[5] C. T. Hill and D. N. Schramm, The Ultrahigh-Energy Cosmic Ray Spectrum , Phys.Rev. D31 
(1985) 564. 

[6] V. Berezinsky, A. Gazizov, and S. Grigorieva, On astrophysical solution to ultrahigh-energy 
cosmic rays, Phys.Rev. D74 (2006) 043005, [hep-ph/0204357]. 

[7] R. Aloisio, V. Berezinsky, P. Blasi, A. Gazizov, S. Grigorieva, et al., A dip in the UHECR 
spectrum and the transition from galactic to extragalactic cosmic rays, Astropart.Phys. 27 
(2007) 76-91, [astro-ph/0608219]. 


- 21 - 







[8] R. Aloisio, V. Berezinsky, and A. Gazizov, Ultra High Energy Cosmic Rays: The disappointing 
model, Astropart.Phys. 34 (2011) 620-626, [arXiv:0907.5194]. 

[9] Pierre Auger Collaboration, A. Aab et al., Depth of maximum of air-shower profiles at the 
Pierre Auger Observatory. II. Composition implications, Phys.Rev. D90 (2014), no. 12 122006, 
[arXiv: 1409.5083], 

[10] A. M. Taylor, M. Aiders, and F. A. Aharonian, Need for a local source of ultrahigh-energy 
cosmic-ray nuclei, Phys.Rev. D84 (2011) 105007, [arXiv: 1107.2055], 

[11] R. Aloisio, V. Berezinsky, and P. Blasi, Ultra high energy cosmic rays: implications of Auger 
data for source spectra and chemical composition, JCAP 1410 (2014), no. 10 020, 

[arXiv: 1312.7459], 

[12] K. Fang, K. Kotera, and A. V. Olinto, Ultrahigh Energy Cosmic Ray Nuclei from Extragalactic 
Pulsars and the effect of their Galactic counterparts, JCAP 1303 (2013) 010, 

[arXiv: 1302.4482], 

[13] A. M. Taylor, UHECR Composition Models, Astropart.Phys. 54 (2014) 48-53, 

[arXiv: 1401.0199], 

[14] C. J. T. Peixoto, V. de Souza, and P. L. Biermann, Cosmic rays: the spectrum and chemical 
composition from 10 10 to 10 20 eV, arXiv: 1502.00305. 

[15] N. Globus, D. Allard, and E. Parizot, A complete model of the CR spectrum and composition 
across the Galactic to Extragalactic transition, arXiv : 1505.01377. 

[16] A. M. Taylor, M. Alders, and D. Hooper, Indications of Negative Evolution for the Sources of 
the Highest Energy Cosmic Rays, arXiv : 1505.06090. 

[17] M. Unger, G. R. Farrar, and L. A. Anchordoqui, Origin of the ankle in the ultra-high energy 
cosmic ray spectrum and of the extragalactic protons below it, arXiv : 1505.02153. 

[18] E. Armengaud, G. Sigl, T. Beau, and F. Miniati, CRPropa: a numerical tool for the 
propagation of UHE cosmic rays, "/-rays and neutrinos, Astropart.Phys. 28 (2007) 463-471, 
[astro-ph/0603675], 

[19] K.-H. Kampert, J. Kulbartz, L. Maccione, N. Nierstenhoefer, P. Schiffer, et al., CRPropa 2.0 
a Public Framework for Propagating High Energy Nuclei, Secondary Gamma Rays and 
Neutrinos, Astropart.Phys. 42 (2013) 41-51, [arXiv: 1206.3132], 

[20] R. Alves Batista, M. Erdmann, C. Evoli, K.-H. Kampert, D. Kuempel, G. Muller, G. Sigl, 

A. van Vliet, D. Walz, and T. Winchen, CRPropa 3, in preparation. 

[21] R. Aloisio, D. Boncioli, A. Grillo, S. Petrera, and F. Salamida, SimProp: a Simulation Code for 
Ultra High Energy Cosmic Ray Propagation, JCAP 1210 (2012) 007, [arXiv: 1204.2970]. 

[22] O. Kalashev and E. Kido, Simulations of Ultra-High-Energy Cosmic Rays propagation , 
J.Exp.Theor.Phys. 120 (2015), no. 5 790-797, [arXiv: 1406.0735], 

[23] M. De Domenico, HERMES: Simulating the Propagation of Ultra-High Energy Cosmic Rays, 
European Physical Journal Plus 128 (Aug., 2013) 99, [arXiv : 1305.4364], 

[24] D. Allard, E. Parizot, E. Khan, S. Goriely, and A. Olinto, UHE nuclei propagation and the 
interpretation of the ankle in the cosmic-ray spectrum, Astron.Astrophys. 443 (2005) L29-L32, 
[astr o-ph/0505566] . 

[25] D. Hooper, S. Sarkar, and A. M. Taylor, The intergalactic propagation of ultra-high energy 
cosmic ray nuclei, Astropart.Phys. 27 (2007) 199-212, [astro-ph/0608085], 

[26] R. Aloisio, V. Berezinsky, and S. Grigorieva, Analytic calculations of the spectra of ultra high 
energy cosmic ray nuclei. II. The general case of background radiation, Astropart.Phys. 41 
(2013) 94-107, [arXiv: 1006.2484], 


22 


[27] R. Gilmore, R. Somerville, J. Primack, and A. Dominguez, Semi-analytic modeling of the EBL 
and consequences for extragalactic gamma-ray spectra , Mon.Not.Roy.Astron.Soc. 422 (2012) 
3189, [arXiv: 1104.0671], 

[28] R. Aloisio and V. Berezinsky, Diffusive propagation of UHECR and the propagation theorem , 
Astrophys.J. 612 (2004) 900-913, [astro-ph/0403095]. 

[29] Pierre Auger Collaboration, P. Abreu et al., Bounds on the density of sources of ultra-high 
energy cosmic rays from the Pierre Auger Observatory , JCAP 1305 (2013), no. 05 009, 

[arXiv: 1305.1576], 

[30] M. Lemoine, Extra-galactic magnetic fields and the second knee in the cosmic-ray spectrum , 
Phys.Rev. D71 (2005) 083007, [astro-ph/0411173]. 

[31] G. Sigl, Non-Universal Spectra of Ultrahigh Energy Cosmic Ray Primaries and Secondaries in 
a Structured Universe, Phys.Rev. D75 (2007) 103001, [astro-ph/0703403]. 

[32] N. Globus, D. Allard, and E. Parizot, Propagation of high-energy cosmic rays in extragalactic 
turbulent magnetic fields: resulting energy spectrum and composition, Astron.Astrophys. 479 
(2008) 97, [arXiv: 0709.1541], 

[33] K. Kotera and M. Lemoine, Inhomogeneous extragalactic magnetic fields and the second knee in 
the cosmic ray spectrum , Phys.Rev. D77 (2008) 023005, [arXiv : 0706.1891], 

[34] S. Mollerach and E. Roulet, Magnetic diffusion effects on the ultra-high energy cosmic ray 
spectrum and composition, JCAP 1310 (2013) 013, [arXiv : 1305.6519]. 

[35] R. Alves Batista and G. Sigl, Diffusion of cosmic rays at EeV energies in inhomogeneous 
extragalactic magnetic fields, JCAP 1411 (2014), no. 11 031, [arXiv : 1407.6150]. 

[36] R. Alves Batista, M. Erdmann, C. Evoli, K.-H. Kampert, D. Kuempel, et al., CRPropa 3.0 - a 
Public Framework for Propagating UHE Cosmic Rays through Galactic and Extragalactic 
Space, arXiv: 1307.2643. 

[37] R. Alves Batista, M. Erdmann, C. Evoli, K.-H. Kampert, D. Kuempel, et al., CRPropa: a 
public framework to propagate UHECRs in the universe, arXiv : 1411.2259. 

[38] R. Alves Batista, M. Erdmann, C. Evoli, K.-H. Kampert, D. Kuempel, et al., Cosmic ray 
propagation with CRPropa 3, J.Phys.Conf.Ser. 608 (2015), no. 1 012076, [arXiv : 1410.5323]. 

[39] G. Blumenthal, Energy loss of high-energy cosmic rays in pair-producing collisions with 
ambient photons, Phys.Rev. D1 (1970) 1596-1602. 

[40] A. Miicke, J. Rachen, R. Engel, R. Protheroe, and T. Stanev, On photohadronic processes in 
astrophysical environments, Publ.Astron.Soc.Austral. 16 (1999) 160, [astro-ph/9808279]. 

[41] A. J. Koning, S. Hilaire, and M. C. Duijvestijn, TALYS: Comprehensive Nuclear Reaction 
Modeling, in International Conference on Nuclear Data for Science and Technology (R. C. 
Haight, M. B. Chadwick, T. Kawano, and P. Talou, eds.), vol. 769 of American Institute of 
Physics Conference Series, pp. 1154-1159, May, 2005. 

[42] J. Rachen, Interaction Processes and Statistical Properties of the Propagation of Cosmic Rays 
in Photon Backgrounds. PhD thesis, University of Bonn, 1996. 

[43] N. Nierstenhoefer, On the Origin and Propagation of Ultra-High Energy Cosmic Rays. PhD 
thesis, University of Wuppertal, 2011. 

[44] M. V. Kossov, Approximation of photonuclear interaction cross-sections, European Physical 
Journal A 14 (2002) 377-392. 

[45] L. A. Kulchitskii, Y. M. Volkov, J. P. Ostriker, and V. I. Ogurtsov, Energy levels of 7 Li 
observed in its photoemission, Izvestiya Rossiiskoi Akademii Nauk. Seriya Fizicheskaya 27 
(1963) 1412. 


- 23 


[46] V. V. Varlamov, Photonuclear data, photodisintegration of lithium, evaluated cross sections of 
channels and reactions , tech, rep., Moscow State University, Moscow, Russia, 1986. 

[47] T. M. Kneiske, T. Bretz, K. Mannheim, and D. Hartmann, Implications of cosmological 
gamma-ray absorption. 2. Modification of gamma-ray spectra, Astron.Astrophys. 413 (2004) 
807-815, [astro-ph/030914l]. 

[48] F. W. Stecker, M. Malkan, and S. Scully, Intergalactic Photon Spectra from the Far IR to the 
UV Lyman Limit for 0 < z < 6 and the Optical Depth of the Universe to High Energy 
Gamma-Rays, Astrophys.J. 648 (2006) 774-783, [astro-ph/0510449]. 

[49] F. W. Stecker, M. Malkan, and S. Scully, Corrected Table for the Parametric Coefficients for 
the Optical Depth of the Universe to Gamma-rays at Various Redshifts, Astrophys.J. 658 
(2007) 1392, [astro-ph/0612048]. 

[50] A. Franceschini, G. Rodighiero, and M. Vaccari, The extragalactic optical-infrared background 
radiations, their time evolution and the cosmic photon-photon opacity, Astron.Astrophys. 487 
(2008) 837, [arXiv:0805.1841], 

[51] J. D. Finke, S. Razzaque, and C. D. Dermer, Modeling the Extragalactic Background Light from 
Stars and Dust, Astrophys.J. 712 (2010) 238-249, [arXiv: 0905.1115]. 

[52] A. Dominguez, J. Primack, D. Rosario, F. Prada, R. Gilmore, et al., Extragalactic Background 
Light Inferred from AEGIS Galaxy SED-type Fractions, Mon.Not.Roy.Astron.Soc. 410 (2011) 
2556, [arXiv: 1007.1459], 

[53] D. Boncioli, Study of Extragalactic Propagation of Cosmic Rays. Applications to Pierre Auger 
Observatory Data. PhD thesis, University of Roma Tor Vergata, 2011. 

[54] J. Puget, F. Stecker, and J. Bredekamp, Photonuclear Interactions of Ultrahigh-Energy Cosmic 
Rays and their Astrophysical Consequences, Astrophys.J. 205 (1976) 638-654. 

[55] F. Stecker and M. Salamon, Photodisintegration of ultrahigh-energy cosmic rays: A New 
determination, Astrophys.J. 512 (1999) 521-526, [astro-ph/9808110]. 

[56] R. Aloisio, V. Berezinsky, and S. Grigorieva, Analytic calculations of the spectra of ultra-high 
energy cosmic ray nuclei. I. The case of CMB radiation, Astropart.Phys. 41 (2013) 73-93, 
[arXiv: 0802.4452], 

[57] R. Aloisio, D. Boncioli, A. Di Matteo, A. Grillo, S. Petrera, et al., Propagation of UHECRs in 
cosmological backgrounds: some results from SimProp, arXiv : 1307.3895. 

[58] R. Aloisio, D. Boncioli, A. di Matteo, A. Grillo, S. Petrera, et al., SimProp v2r2: a Monte 
Carlo simulation to compute cosmogenic neutrino fluxes, arXiv: 1505.01347. 

[59] Pierre Auger Collaboration, A. Aab et al., Depth of maximum of air-shower profiles at the 
Pierre Auger Observatory. I. Measurements at energies above 10 ll 8 eV, Phys.Rev. D90 (2014), 
no. 12 122005, [arXiv: 1409.4809], 

[60] Pierre Auger Collaboration, A. Aab et al., The Pierre Auger Observatory: Contributions to 
the 33rd International Cosmic Ray Conference (ICRC 2013), arXiv : 1307.5059. 

[61] T. Pierog, I. Karpenko, J. Katzy, E. Yatsenko, and K. Werner, EPOS LHC : test of collective 
hadronization with LHC data, arXiv : 1306.0121. 

[62] D. Allard, Extragalactic propagation of ultrahigh energy cosmic-rays, Astropart.Phys. 39-40 
(2012) 33-43, [arXiv: 1111.3290], 

[63] E. Dwek and F. Krennrich, The Extragalactic Background Light and the Gamma-ray Opacity of 
the Universe, Astropart.Phys. 43 (2013) 112-133, [arXiv : 1209.4661]. 

[64] R. Aloisio, D. Boncioli, A. di Matteo, A. Grillo, S. Petrera, et al., Cosmogenic neutrinos and 
ultra-high energy cosmic ray models, arXiv : 1505.04020. 


- 24 - 


[65] M. Malkan and F. Stecker, An empirically based model for predicting infrared luminosity 
functions, deep infrared galaxy counts and the diffuse infrared background, Astrophys.J. 555 
(2001) 641-649, [astro-ph/0009500]. 

[66] HEGRA Collaboration, F. Aharonian et al., Is the giant radio galaxy M87 a TeV gamma-ray 
emitter?, Astron.Astrophys. 403 (2003) L1-L6, [astro-ph/0302155]. 

[67] A. Franceschini, H. Aussel, C. Cesarsky, D. Elbaz, and D. Fadda, A long-wavelength view on 
galaxy evolution from deep surveys by the infrared space observatory, Astron.Astrophys. 378 
(2001) 1-29, [astro-ph/0108292], 

[68] E. Khan et al., Photodisintegration of ultra-high-energy cosmic rays revisited, Astroparticle 
Physics 23 (Mar., 2005) 191-201, [astro-ph/0412109]. 

[69] M. Chadwick et al., Handbook on photonuclear data for applications: cross-sections and 
spectra, tech, rep., IAEA, 2000. 

[70] T. Belgya, O. Bersillon, R. C. Noy, T. Fukahori, G. Zhigang, S. Goriely, M. Herman, 

A. Ignatyuk, S. Kailas, A. Koning, P. Oblozinsky, V. Plujko, and P. Young, Handbook for 
calculations of nuclear reaction data, RIPL-2, Tech. Rep. 1506, IAEA, 2006. 

[71] A. J. Koning and S. Goriely. private communication, 2015. 

[72] A. Koning, S. Hilaire, and S. Goriely, TALYS 1.6 User Manual. 

[73] B. Ishkhanov et al., Cross sections of photon absorption by nuclei with nucleon numbers 12 - 
65, Rept: Moscow State Univ. Inst, of Nucl. Phys. Reports 2002 (2002) 27-271. 

[74] S. Afanas’ev and A. Khodyachikh, On the mechanism of formation of excited states of the 8 Be 
nucleus in the reaction 12 (7(7, 3 a), Physics of Atomic Nuclei 71 (2008), no. 11 1827-1838. 


- 25 - 


